Code
library(ggplot2)
library(pander)
library(emmeans)
library(tidyverse)
library(here)
library(ggpubr)
library(ggridges)
library(readxl)
library(magrittr)
library(patchwork)
library(here)library(ggplot2)
library(pander)
library(emmeans)
library(tidyverse)
library(here)
library(ggpubr)
library(ggridges)
library(readxl)
library(magrittr)
library(patchwork)
library(here)Analysis preps (colours, options, auxiliary functions)
cbb_palette <- c(
"#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442",
"#0072B2", "#D55E00", "#CC79A7"
)
## Avoid table wrapping
options(width = 800)
set.seed(999)
source(here("PhillipsKarpEtal_data/funs.r"))Below results repeat simulations from the paper of Phillips et al. (2018) PLoS Biol 21(6): e3002129. Original results simplified by removing the post-hoc and the (pooled) t-test methods of analysing data. Generation of simulated data was modified to use a formal linear model defined as:
\[\mathbf{y} = \mathbf{X}\mathbf{b} + \mathbf{e}\]
with \(\mathbf{b}=(\mu,\beta_{Treat},\beta_{Sex},\beta_{Treat:Sex})\) and \(\mathbf{e} \sim \mathcal{N}(\mathbf{0}, \mathbf{s})\), with \(\mathbf{s}=\sigma\sqrt{\mathbf{1}+\mathbf{X}_{Sex}\cdot\alpha}\) (\(\alpha\) quantifies the amount of heteroscedasticity in sex variances, with SD \(\sigma_{m}=\sigma_{f}\sqrt{1+\alpha}\), or male variance being \(\alpha\cdot100\%\) greater than female). Regression coefficients (\(\mathbf{b}\)) are compatible with contrasts (sex/treatment effects) used in Phillips et al. (2023).
Assuming homoscedasticity of sex variances
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+0
all_effs_sex_none <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_none <- all_effs_sex_none %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("none")) %>%
mutate(heterosc = rep("none"))
calc_power_sex_none %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1
all_effs_sex_small <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_small <- all_effs_sex_small %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("small")) %>%
mutate(heterosc = rep("none"))
calc_power_sex_small %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1
all_effs_sex_large <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_large <- all_effs_sex_large %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(heterosc = rep("none"))
calc_power_sex_large %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)Plots combining all scenarios from Fig. 1
all_sex_effs_sex_only <- rbind(calc_power_sex_none, calc_power_sex_small,
calc_power_sex_large)
all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))
all_sex_effs_sex_only %>%
filter(Effect %in% c("Treatment")) %>%
mutate(sex_eff_size = factor(sex_eff_size,
levels = c("none", "small", "large"))) %>%
mutate(method = recode(method,
"ANOVA" = "Factorial")) %>%
ggplot(aes(x = treat_es2, y = p_value, colour = sex_eff_size,
group = interaction(sex_eff_size, method))) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(linetype = method)) +
theme_classic() +
ylab("Statistical power") +
xlab("Main treatment effect") +
scale_colour_manual(name = "Sex effect size", values = cbb_palette) +
scale_linetype_discrete(name = "Stat method") + ylim(0, 1)all_sex_effs_sex_only %>%
# filter(!Effect %in% c("T-Test Treatment Eff")) %>%
mutate(sex_eff_size = factor(sex_eff_size,
levels = c("none", "small", "large"))) %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect,
group = interaction(sex_eff_size, Effect))) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(linetype = sex_eff_size)) +
theme_classic() +
ylab("Statistical power") +
xlab("Main treatment effect") +
scale_colour_manual(name = "Effect", values = cbb_palette) +
scale_linetype_discrete(name = "Sex effect size") + ylim(0, 1)We assume several scenarios with \(\alpha\) equal to 0 (the above scenarios), 1 and 2, which relates to 100% or 200% increase in male variance, in relation to female variance (or ~44% and ~73% increase in male SD). The direction of sex bias is arbitrary.
Low heteroscedasticity
alpha <- 1no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha
all_effs_sex_none_lowh <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_none_lowh <- all_effs_sex_none_lowh %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("none")) %>%
mutate(heterosc = rep("low"))
calc_power_sex_none_lowh %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha
all_effs_sex_small_lowh <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_small_lowh <- all_effs_sex_small_lowh %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("small")) %>%
mutate(heterosc = rep("low"))
calc_power_sex_small_lowh %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha
all_effs_sex_large_lowh <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_large_lowh <- all_effs_sex_large_lowh %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(heterosc = rep("low"))
calc_power_sex_large_lowh %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)Large heteroscedasticity
alpha <- 2no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha
all_effs_sex_none_largeh <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_none_largeh <- all_effs_sex_none_largeh %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("none")) %>%
mutate(heterosc = rep("large"))
calc_power_sex_none_largeh %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha
all_effs_sex_small_largeh <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_small_largeh <- all_effs_sex_small_largeh %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("small")) %>%
mutate(heterosc = rep("large"))
calc_power_sex_small_largeh %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha
all_effs_sex_large_largeh <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_large_largeh <- all_effs_sex_large_largeh %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(heterosc = rep("large"))
calc_power_sex_large_largeh %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)It should be of note that evolutionary biology has seen heteroscedasticities exceeding 100% in variance ration between sexes. Combined results below show the impact of heteroscedasticity.
all_sex_effs_sex_only <- rbind(calc_power_sex_small,
calc_power_sex_large,
calc_power_sex_small_lowh,
calc_power_sex_large_lowh,
calc_power_sex_small_largeh,
calc_power_sex_large_largeh)
all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))
all_sex_effs_sex_only %>%
filter(Effect %in% c("Treatment")) %>%
mutate(sex_eff_size = factor(sex_eff_size,
levels = c("small", "large"))) %>%
mutate(method = recode(method,
"ANOVA" = "Factorial")) %>%
mutate(heterosc = factor(heterosc,
levels = c("none", "low", "large"))) %>%
ggplot(aes(x = treat_es2, y = p_value, colour = sex_eff_size,
group = interaction(sex_eff_size, heterosc))) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(linetype = heterosc)) +
theme_classic() +
ylab("Statistical power") +
xlab("Main treatment effect") +
scale_colour_manual(name = "Sex effect size", values = cbb_palette) +
scale_linetype_discrete(name = "Heteroscedasticity") + ylim(0, 1)plotA <- all_sex_effs_sex_only %>%
filter(sex_eff_size %in% c("large")) %>%
# mutate(sex_eff_size = factor(sex_eff_size,
# levels = c("none", "small", "large"))) %>%
mutate(heterosc = factor(heterosc,
levels = c("none", "low", "large"))) %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(linetype = heterosc)) +
theme_classic() +
ylab("Statistical power") +
xlab("Main treatment effect") +
scale_colour_manual(name = "Effect", values = cbb_palette) +
scale_linetype_discrete(name = "Heteroscedasticity") + ylim(0, 1)
plotABelow simulations were only slightly modified. We assumed a moderate (+0.5) overall effect of sex and a varying interaction increasing the magnitude of sex difference in the treated group (+0, +0.5, +1).
Assuming homoscedasticity of sex variances
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+0
all_effs_sex_ix_none <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_none <- all_effs_sex_ix_none %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("small")) %>%
mutate(heterosc = rep("none")) %>%
mutate(ix_eff_size = rep("none"))
calc_power_sex_ix_none %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0.5
sex_sd <- 1
all_effs_sex_ix_small <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_small <- all_effs_sex_ix_small %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("small")) %>%
mutate(heterosc = rep("none")) %>%
mutate(ix_eff_size = rep("small"))
calc_power_sex_ix_small %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 1
sex_sd <- 1
all_effs_sex_ix_large <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_large <- all_effs_sex_ix_large %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(heterosc = rep("none")) %>%
mutate(ix_eff_size = rep("large"))
calc_power_sex_ix_large %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)Plots combining all scenarios from Fig. 1
all_sex_effs_sex_only <- rbind(calc_power_sex_ix_none, calc_power_sex_ix_small,
calc_power_sex_ix_large)
all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))
all_sex_effs_sex_only %>%
filter(Effect %in% c("Treatment")) %>%
mutate(ix_eff_size = factor(ix_eff_size,
levels = c("none", "small", "large"))) %>%
mutate(method = recode(method,
"ANOVA" = "Factorial")) %>%
ggplot(aes(x = treat_es2, y = p_value, colour = ix_eff_size,
group = interaction(ix_eff_size, method))) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(linetype = method)) +
theme_classic() +
ylab("Statistical power") +
xlab("Main treatment effect") +
scale_colour_manual(name = "Ix effect size", values = cbb_palette) +
scale_linetype_discrete(name = "Stat method") + ylim(0, 1)all_sex_effs_sex_only %>%
# filter(!Effect %in% c("T-Test Treatment Eff")) %>%
mutate(ix_eff_size = factor(ix_eff_size,
levels = c("none", "small", "large"))) %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect,
group = interaction(ix_eff_size, Effect))) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(linetype = ix_eff_size)) +
theme_classic() +
ylab("Statistical power") +
xlab("Main treatment effect") +
scale_colour_manual(name = "Effect", values = cbb_palette) +
scale_linetype_discrete(name = "Ix effect size") + ylim(0, 1)We assume several scenarios with \(\alpha\) equal to 0 (the above scenarios), 1 and 2, which relates to 100% or 200% increase in male variance, in relation to female variance (or ~44% and ~73% increase in male SD). The direction of sex bias is arbitrary.
Low heteroscedasticity
alpha <- 1no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha
all_effs_sex_ix_none_lowh <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_none_lowh <- all_effs_sex_ix_none_lowh %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("small")) %>%
mutate(heterosc = rep("low")) %>%
mutate(ix_eff_size = rep("none"))
calc_power_sex_ix_none_lowh %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0.5
sex_sd <- 1+alpha
all_effs_sex_ix_small_lowh <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_small_lowh <- all_effs_sex_ix_small_lowh %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("small")) %>%
mutate(heterosc = rep("low")) %>%
mutate(ix_eff_size = rep("small"))
calc_power_sex_ix_small_lowh %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 1
sex_sd <- 1+alpha
all_effs_sex_ix_large_lowh <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_large_lowh <- all_effs_sex_ix_large_lowh %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("small")) %>%
mutate(heterosc = rep("low")) %>%
mutate(ix_eff_size = rep("large"))
calc_power_sex_ix_large_lowh %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)Large heteroscedasticity
alpha <- 2no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha
all_effs_sex_ix_none_largeh <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_none_largeh <- all_effs_sex_ix_none_largeh %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("small")) %>%
mutate(heterosc = rep("large")) %>%
mutate(ix_eff_size = rep("none"))
calc_power_sex_ix_none_largeh %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0.5
sex_sd <- 1+alpha
all_effs_sex_ix_small_largeh <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_small_largeh <- all_effs_sex_ix_small_largeh %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("small")) %>%
mutate(heterosc = rep("large")) %>%
mutate(ix_eff_size = rep("small"))
calc_power_sex_ix_small_largeh %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 0.5
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 1
sex_sd <- 1+alpha
all_effs_sex_ix_large_largeh <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_large_largeh <- all_effs_sex_ix_large_largeh %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(heterosc = rep("large")) %>%
mutate(ix_eff_size = rep("large"))
calc_power_sex_ix_large_largeh %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)Power consequences of heteroscedasticities with varying (magnitude-only) interaction.
all_sex_effs_sex_only <- rbind(calc_power_sex_ix_small,
calc_power_sex_ix_large,
calc_power_sex_ix_small_lowh,
calc_power_sex_ix_large_lowh,
calc_power_sex_ix_small_largeh,
calc_power_sex_ix_large_largeh)
all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))
all_sex_effs_sex_only %>%
filter(Effect %in% c("Treatment")) %>%
mutate(ix_eff_size = factor(ix_eff_size,
levels = c("small", "large"))) %>%
mutate(method = recode(method,
"ANOVA" = "Factorial")) %>%
mutate(heterosc = factor(heterosc,
levels = c("none", "low", "large"))) %>%
ggplot(aes(x = treat_es2, y = p_value, colour = ix_eff_size,
group = interaction(ix_eff_size, heterosc))) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(linetype = heterosc)) +
theme_classic() +
ylab("Statistical power") +
xlab("Main treatment effect") +
scale_colour_manual(name = "Ix effect size", values = cbb_palette) +
scale_linetype_discrete(name = "Heteroscedasticity") + ylim(0, 1)plotB <- all_sex_effs_sex_only %>%
filter(ix_eff_size %in% c("large")) %>%
# mutate(sex_eff_size = factor(sex_eff_size,
# levels = c("none", "small", "large"))) %>%
mutate(heterosc = factor(heterosc,
levels = c("none", "low", "large"))) %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(linetype = heterosc)) +
theme_classic() +
ylab("Statistical power") +
xlab("Main treatment effect") +
scale_colour_manual(name = "Effect", values = cbb_palette) +
scale_linetype_discrete(name = "Heteroscedasticity") + ylim(0, 1)
plotBBelow simulations were only slightly modified. We assumed a strong (-1) overall effect of sex and a varying interaction resulting of crossing of sex effects on the treatment effect (0, -1, -2).
Assuming homoscedasticity of sex variances
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+0
all_effs_sex_ix_none <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_none <- all_effs_sex_ix_none %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("small")) %>%
mutate(heterosc = rep("none")) %>%
mutate(ix_eff_size = rep("none"))
calc_power_sex_ix_none %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -1
sex_sd <- 1
all_effs_sex_ix_small <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_small <- all_effs_sex_ix_small %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("small")) %>%
mutate(heterosc = rep("none")) %>%
mutate(ix_eff_size = rep("small"))
calc_power_sex_ix_small %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -2
sex_sd <- 1
all_effs_sex_ix_large <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_large <- all_effs_sex_ix_large %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(heterosc = rep("none")) %>%
mutate(ix_eff_size = rep("large"))
calc_power_sex_ix_large %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)Plots combining all scenarios from Fig. 1
all_sex_effs_sex_only <- rbind(calc_power_sex_ix_none, calc_power_sex_ix_small,
calc_power_sex_ix_large)
all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))
all_sex_effs_sex_only %>%
filter(Effect %in% c("Treatment")) %>%
mutate(ix_eff_size = factor(ix_eff_size,
levels = c("none", "small", "large"))) %>%
mutate(method = recode(method,
"ANOVA" = "Factorial")) %>%
ggplot(aes(x = treat_es2, y = p_value, colour = ix_eff_size,
group = interaction(ix_eff_size, method))) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(linetype = method)) +
theme_classic() +
ylab("Statistical power") +
xlab("Main treatment effect") +
scale_colour_manual(name = "Ix effect size", values = cbb_palette) +
scale_linetype_discrete(name = "Stat method") + ylim(0, 1)all_sex_effs_sex_only %>%
# filter(!Effect %in% c("T-Test Treatment Eff")) %>%
mutate(ix_eff_size = factor(ix_eff_size,
levels = c("none", "small", "large"))) %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect,
group = interaction(ix_eff_size, Effect))) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(linetype = ix_eff_size)) +
theme_classic() +
ylab("Statistical power") +
xlab("Main treatment effect") +
scale_colour_manual(name = "Effect", values = cbb_palette) +
scale_linetype_discrete(name = "Ix effect size") + ylim(0, 1)We assume several scenarios with \(\alpha\) equal to 0 (the above scenarios), 1 and 2, which relates to 100% or 200% increase in male variance, in relation to female variance (or ~44% and ~73% increase in male SD). The direction of sex bias is arbitrary.
Low heteroscedasticity
alpha <- 1no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha
all_effs_sex_ix_none_lowh <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_none_lowh <- all_effs_sex_ix_none_lowh %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("small")) %>%
mutate(heterosc = rep("low")) %>%
mutate(ix_eff_size = rep("none"))
calc_power_sex_ix_none_lowh %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -1
sex_sd <- 1+alpha
all_effs_sex_ix_small_lowh <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_small_lowh <- all_effs_sex_ix_small_lowh %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("small")) %>%
mutate(heterosc = rep("low")) %>%
mutate(ix_eff_size = rep("small"))
calc_power_sex_ix_small_lowh %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -2
sex_sd <- 1+alpha
all_effs_sex_ix_large_lowh <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_large_lowh <- all_effs_sex_ix_large_lowh %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("small")) %>%
mutate(heterosc = rep("low")) %>%
mutate(ix_eff_size = rep("large"))
calc_power_sex_ix_large_lowh %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)Large heteroscedasticity
alpha <- 2no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha
all_effs_sex_ix_none_largeh <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_none_largeh <- all_effs_sex_ix_none_largeh %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("small")) %>%
mutate(heterosc = rep("large")) %>%
mutate(ix_eff_size = rep("none"))
calc_power_sex_ix_none_largeh %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -1
sex_sd <- 1+alpha
all_effs_sex_ix_small_largeh <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_small_largeh <- all_effs_sex_ix_small_largeh %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("small")) %>%
mutate(heterosc = rep("large")) %>%
mutate(ix_eff_size = rep("small"))
calc_power_sex_ix_small_largeh %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -2
sex_sd <- 1+alpha
all_effs_sex_ix_large_largeh <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_large_largeh <- all_effs_sex_ix_large_largeh %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, ix_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(heterosc = rep("large")) %>%
mutate(ix_eff_size = rep("large"))
calc_power_sex_ix_large_largeh %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)Power consequences of heteroscedasticities with varying (direction-only) interaction. Power anomalies indicate pretty low inferential abilities of such models at low sample sizes.
all_sex_effs_sex_only <- rbind(calc_power_sex_ix_small,
calc_power_sex_ix_large,
calc_power_sex_ix_small_lowh,
calc_power_sex_ix_large_lowh,
calc_power_sex_ix_small_largeh,
calc_power_sex_ix_large_largeh)
all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))
all_sex_effs_sex_only %>%
filter(Effect %in% c("Treatment")) %>%
mutate(ix_eff_size = factor(ix_eff_size,
levels = c("small", "large"))) %>%
mutate(method = recode(method,
"ANOVA" = "Factorial")) %>%
mutate(heterosc = factor(heterosc,
levels = c("none", "low", "large"))) %>%
ggplot(aes(x = treat_es2, y = p_value, colour = ix_eff_size,
group = interaction(ix_eff_size, heterosc))) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(linetype = heterosc)) +
theme_classic() +
ylab("Statistical power") +
xlab("Main treatment effect") +
scale_colour_manual(name = "Ix effect size", values = cbb_palette) +
scale_linetype_discrete(name = "Heteroscedasticity") + ylim(0, 1)all_sex_effs_sex_only %>%
filter(ix_eff_size %in% c("large")) %>%
# mutate(sex_eff_size = factor(sex_eff_size,
# levels = c("none", "small", "large"))) %>%
mutate(heterosc = factor(heterosc,
levels = c("none", "low", "large"))) %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(linetype = heterosc)) +
theme_classic() +
ylab("Statistical power") +
xlab("Main treatment effect") +
scale_colour_manual(name = "Effect", values = cbb_palette) +
scale_linetype_discrete(name = "Heteroscedasticity") + ylim(0, 1)Below simulation code generates power estimates for the following sample size scenarios: 5M+5F in each treatment group (original scenario, totql sample 20 individuals), 10M+10F (N = 40), 50M+50F (N = 200).
First - low heteroscedasticity:
alpha <- 1
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha
all_effs_sex_large_lowh_5 <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_large_lowh_5 <- all_effs_sex_large_lowh_5 %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(ix_eff_size = rep("none")) %>%
mutate(sample_size = rep("5M+5F")) %>%
mutate(heterosc = rep("low"))
calc_power_sex_large_lowh_5 %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)no_per_gp <- 10
all_effs_sex_large_lowh_10 <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_large_lowh_10 <- all_effs_sex_large_lowh_10 %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(ix_eff_size = rep("none")) %>%
mutate(sample_size = rep("10M+10F")) %>%
mutate(heterosc = rep("low"))
calc_power_sex_large_lowh_10 %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)no_per_gp <- 50
all_effs_sex_large_lowh_50 <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_large_lowh_50 <- all_effs_sex_large_lowh_50 %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(ix_eff_size = rep("none")) %>%
mutate(sample_size = rep("50M+50F")) %>%
mutate(heterosc = rep("low"))
calc_power_sex_large_lowh_50 %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)Large heteroscedasticity:
alpha <- 2
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha
all_effs_sex_large_largeh_5 <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_large_largeh_5 <- all_effs_sex_large_largeh_5 %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(ix_eff_size = rep("none")) %>%
mutate(sample_size = rep("5M+5F")) %>%
mutate(heterosc = rep("large"))
calc_power_sex_large_largeh_5 %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)no_per_gp <- 10
all_effs_sex_large_largeh_10 <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_large_largeh_10 <- all_effs_sex_large_largeh_10 %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(ix_eff_size = rep("none")) %>%
mutate(sample_size = rep("10M+10F")) %>%
mutate(heterosc = rep("large"))
calc_power_sex_large_largeh_10 %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)no_per_gp <- 50
all_effs_sex_large_largeh_50 <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_large_largeh_50 <- all_effs_sex_large_largeh_50 %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(ix_eff_size = rep("none")) %>%
mutate(sample_size = rep("50M+50F")) %>%
mutate(heterosc = rep("large"))
calc_power_sex_large_largeh_50 %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)Identical sample size scenarios.
First - low heteroscedasticity
alpha <- 1
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -2
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -2
sex_sd <- 1+alpha
all_effs_sex_ix_large_lowh_5 <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_large_lowh_5 <- all_effs_sex_ix_large_lowh_5 %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(ix_eff_size = rep("large")) %>%
mutate(sample_size = rep("5M+5F")) %>%
mutate(heterosc = rep("low"))
calc_power_sex_ix_large_lowh_5 %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)no_per_gp <- 10
all_effs_sex_ix_large_lowh_10 <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_large_lowh_10 <- all_effs_sex_ix_large_lowh_10 %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(ix_eff_size = rep("large")) %>%
mutate(sample_size = rep("10M+10F")) %>%
mutate(heterosc = rep("low"))
calc_power_sex_ix_large_lowh_10 %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)no_per_gp <- 50
all_effs_sex_ix_large_lowh_50 <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_large_lowh_50 <- all_effs_sex_ix_large_lowh_50 %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(ix_eff_size = rep("large")) %>%
mutate(sample_size = rep("50M+50F")) %>%
mutate(heterosc = rep("low"))
calc_power_sex_ix_large_lowh_50 %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)Large heteroscedasticity
alpha <- 2
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- -2
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -2
sex_sd <- 1+alpha
all_effs_sex_ix_large_largeh_5 <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_large_largeh_5 <- all_effs_sex_ix_large_largeh_5 %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(ix_eff_size = rep("large")) %>%
mutate(sample_size = rep("5M+5F")) %>%
mutate(heterosc = rep("large"))
calc_power_sex_ix_large_largeh_5 %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)no_per_gp <- 10
all_effs_sex_ix_large_largeh_10 <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_large_largeh_10 <- all_effs_sex_ix_large_largeh_10 %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(ix_eff_size = rep("large")) %>%
mutate(sample_size = rep("10M+10F")) %>%
mutate(heterosc = rep("large"))
calc_power_sex_ix_large_largeh_10 %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)no_per_gp <- 50
all_effs_sex_ix_large_largeh_50 <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd
)
)
calc_power_sex_ix_large_largeh_50 <- all_effs_sex_ix_large_largeh_50 %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(ix_eff_size = rep("large")) %>%
mutate(sample_size = rep("50M+50F")) %>%
mutate(heterosc = rep("large"))
calc_power_sex_ix_large_largeh_50 %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(group = Effect)) +
theme_classic() +
ylab("Power") + ylim(0, 1)No interaction, low (left) and large (right) heteroscedasticity:
all_sex_effs_sex_only <- rbind(
calc_power_sex_large_lowh_5,
calc_power_sex_large_lowh_10,
calc_power_sex_large_lowh_50,
calc_power_sex_large_largeh_5,
calc_power_sex_large_largeh_10,
calc_power_sex_large_largeh_50
)
all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))
# plot1 <- all_sex_effs_sex_only %>%
# filter(Effect %in% c("Treatment")) %>%
# filter(heterosc %in% c("low")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(method = recode(method,
# "ANOVA" = "Factorial")) %>%
# ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
# geom_line() +
# theme_classic() +
# ylab("Statistical power") +
# xlab("Main treatment effect") +
# scale_colour_manual(name = "Sample size", values = cbb_palette) +
# # scale_linetype_discrete(name = "Heteroscedasticity") +
# ylim(0, 1)
# plot2 <- all_sex_effs_sex_only %>%
# filter(Effect %in% c("Treatment")) %>%
# filter(heterosc %in% c("large")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(method = recode(method,
# "ANOVA" = "Factorial")) %>%
# ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
# geom_line() +
# theme_classic() +
# ylab("Statistical power") +
# xlab("Main treatment effect") +
# scale_colour_manual(name = "Sample size", values = cbb_palette) +
# # scale_linetype_discrete(name = "Heteroscedasticity") +
# ylim(0, 1)
# ggarrange(plot1, plot2, ncol = 2, nrow = 1, common.legend = TRUE, legend = "right")
plot1 <- all_sex_effs_sex_only %>%
filter(heterosc %in% c("low")) %>%
mutate(sample_size = factor(sample_size,
levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(sex_eff_size = factor(sex_eff_size,
# levels = c("none", "small", "large"))) %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(linetype = sample_size)) +
theme_classic() +
ylab("Statistical power") +
xlab("Main treatment effect") +
scale_colour_manual(name = "Effect", values = cbb_palette) +
scale_linetype_discrete(name = "Sample size") + ylim(0, 1)
plot2 <- all_sex_effs_sex_only %>%
filter(heterosc %in% c("large")) %>%
mutate(sample_size = factor(sample_size,
levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(sex_eff_size = factor(sex_eff_size,
# levels = c("none", "small", "large"))) %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(linetype = sample_size)) +
theme_classic() +
ylab("Statistical power") +
xlab("Main treatment effect") +
scale_colour_manual(name = "Effect", values = cbb_palette) +
scale_linetype_discrete(name = "Sample size") + ylim(0, 1)
ggarrange(plot1, plot2, ncol = 2, nrow = 1, common.legend = TRUE, legend = "right")Interaction of effect direction present, low (left) and high (right) heteroscedasticity:
all_sex_effs_sex_only <- rbind(
calc_power_sex_ix_large_lowh_5,
calc_power_sex_ix_large_lowh_10,
calc_power_sex_ix_large_lowh_50,
calc_power_sex_ix_large_largeh_5,
calc_power_sex_ix_large_largeh_10,
calc_power_sex_ix_large_largeh_50
)
all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))
# plot1 <- all_sex_effs_sex_only %>%
# filter(Effect %in% c("Treatment")) %>%
# filter(heterosc %in% c("low")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(method = recode(method,
# "ANOVA" = "Factorial")) %>%
# ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
# geom_line() +
# theme_classic() +
# ylab("Statistical power") +
# xlab("Main treatment effect") +
# scale_colour_manual(name = "Sample size", values = cbb_palette) +
# # scale_linetype_discrete(name = "Heteroscedasticity") +
# ylim(0, 1)
# plot2 <- all_sex_effs_sex_only %>%
# filter(Effect %in% c("Treatment")) %>%
# filter(heterosc %in% c("large")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(method = recode(method,
# "ANOVA" = "Factorial")) %>%
# ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
# geom_line() +
# theme_classic() +
# ylab("Statistical power") +
# xlab("Main treatment effect") +
# scale_colour_manual(name = "Sample size", values = cbb_palette) +
# # scale_linetype_discrete(name = "Heteroscedasticity") +
# ylim(0, 1)
# ggarrange(plot1, plot2, ncol = 2, nrow = 1, common.legend = TRUE, legend = "right")
plot1 <- all_sex_effs_sex_only %>%
filter(heterosc %in% c("low")) %>%
mutate(sample_size = factor(sample_size,
levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(sex_eff_size = factor(sex_eff_size,
# levels = c("none", "small", "large"))) %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(linetype = sample_size)) +
theme_classic() +
ylab("Statistical power") +
xlab("Main treatment effect") +
scale_colour_manual(name = "Effect", values = cbb_palette) +
scale_linetype_discrete(name = "Sample size") + ylim(0, 1)
plot2 <- all_sex_effs_sex_only %>%
filter(heterosc %in% c("large")) %>%
mutate(sample_size = factor(sample_size,
levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(sex_eff_size = factor(sex_eff_size,
# levels = c("none", "small", "large"))) %>%
ggplot(aes(x = treat_es2, y = p_value, colour = Effect)) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line(aes(linetype = sample_size)) +
theme_classic() +
ylab("Statistical power") +
xlab("Main treatment effect") +
scale_colour_manual(name = "Effect", values = cbb_palette) +
scale_linetype_discrete(name = "Sample size") + ylim(0, 1)
ggarrange(plot1, plot2, ncol = 2, nrow = 1, common.legend = TRUE, legend = "right")Scenario 1: strong sex effect, none to large heteroscedasticity. Assumed power 80% for treatment effect.
alpha <- 0
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha
all_effs_sex_large_noh <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
)
)
calc_power_sex_large_noh <- all_effs_sex_large_noh %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(ix_eff_size = rep("none")) %>%
mutate(heterosc = rep("none"))
alpha <- 1
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha
all_effs_sex_large_lowh <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
)
)
calc_power_sex_large_lowh <- all_effs_sex_large_lowh %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(ix_eff_size = rep("none")) %>%
mutate(heterosc = rep("low"))
alpha <- 2
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0
sex_sd <- 1+alpha
all_effs_sex_large_largeh <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
)
)
calc_power_sex_large_largeh <- all_effs_sex_large_largeh %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(ix_eff_size = rep("none")) %>%
mutate(heterosc = rep("large"))
all_sex_effs_sex_only <- rbind(
calc_power_sex_large_noh,
calc_power_sex_large_lowh,
calc_power_sex_large_largeh
)
all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))
# plot1 <- all_sex_effs_sex_only %>%
# filter(Effect %in% c("Treatment")) %>%
# filter(heterosc %in% c("low")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(method = recode(method,
# "ANOVA" = "Factorial")) %>%
# ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
# geom_line() +
# theme_classic() +
# ylab("Statistical power") +
# xlab("Main treatment effect") +
# scale_colour_manual(name = "Sample size", values = cbb_palette) +
# # scale_linetype_discrete(name = "Heteroscedasticity") +
# ylim(0, 1)
# plot2 <- all_sex_effs_sex_only %>%
# filter(Effect %in% c("Treatment")) %>%
# filter(heterosc %in% c("large")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(method = recode(method,
# "ANOVA" = "Factorial")) %>%
# ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
# geom_line() +
# theme_classic() +
# ylab("Statistical power") +
# xlab("Main treatment effect") +
# scale_colour_manual(name = "Sample size", values = cbb_palette) +
# # scale_linetype_discrete(name = "Heteroscedasticity") +
# ylim(0, 1)
# ggarrange(plot1, plot2, ncol = 2, nrow = 1, common.legend = TRUE, legend = "right")
plot1 <- all_sex_effs_sex_only %>%
# filter(heterosc %in% c("low")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(sex_eff_size = factor(sex_eff_size,
# levels = c("none", "small", "large"))) %>%
ggplot(aes(x = treat_es2, y = ssize, colour = heterosc)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line() +
theme_bw() +
ylab("Sample size") +
xlab("Main treatment effect") +
scale_colour_manual(name = "Heteroscedasticity", values = cbb_palette) + coord_trans(y = "log10") +
theme(text = element_text(size=16))
# plot1
plot1.1 <- all_sex_effs_sex_only %>%
# filter(heterosc %in% c("low")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(sex_eff_size = factor(sex_eff_size,
# levels = c("none", "small", "large"))) %>%
ggplot(aes(x = treat_es2, y = ssize, colour = heterosc)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line() +
theme_bw() +
ylab("Sample size") +
xlab("Main treatment effect") +
scale_colour_manual(name = "Heteroscedasticity", values = cbb_palette) +
xlim(0.3, 0.7) + ylim(5, 70)
# plot1.1Scenario 2: same but interaction of magnitude with more variable sex responding more strongly to the treatment.
alpha <- 0
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0.5
sex_sd <- 1+alpha
all_effs_sex_large_noh <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
)
)
calc_power_sex_large_noh <- all_effs_sex_large_noh %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(ix_eff_size = rep("none")) %>%
mutate(heterosc = rep("none"))
alpha <- 1
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0.5
sex_sd <- 1+alpha
all_effs_sex_large_lowh <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
)
)
calc_power_sex_large_lowh <- all_effs_sex_large_lowh %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(ix_eff_size = rep("none")) %>%
mutate(heterosc = rep("low"))
alpha <- 2
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- 0.5
sex_sd <- 1+alpha
all_effs_sex_large_largeh <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
)
)
calc_power_sex_large_largeh <- all_effs_sex_large_largeh %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(ix_eff_size = rep("none")) %>%
mutate(heterosc = rep("large"))
all_sex_effs_sex_only <- rbind(
calc_power_sex_large_noh,
calc_power_sex_large_lowh,
calc_power_sex_large_largeh
)
all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))
# plot1 <- all_sex_effs_sex_only %>%
# filter(Effect %in% c("Treatment")) %>%
# filter(heterosc %in% c("low")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(method = recode(method,
# "ANOVA" = "Factorial")) %>%
# ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
# geom_line() +
# theme_classic() +
# ylab("Statistical power") +
# xlab("Main treatment effect") +
# scale_colour_manual(name = "Sample size", values = cbb_palette) +
# # scale_linetype_discrete(name = "Heteroscedasticity") +
# ylim(0, 1)
# plot2 <- all_sex_effs_sex_only %>%
# filter(Effect %in% c("Treatment")) %>%
# filter(heterosc %in% c("large")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(method = recode(method,
# "ANOVA" = "Factorial")) %>%
# ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
# geom_line() +
# theme_classic() +
# ylab("Statistical power") +
# xlab("Main treatment effect") +
# scale_colour_manual(name = "Sample size", values = cbb_palette) +
# # scale_linetype_discrete(name = "Heteroscedasticity") +
# ylim(0, 1)
# ggarrange(plot1, plot2, ncol = 2, nrow = 1, common.legend = TRUE, legend = "right")
plot2 <- all_sex_effs_sex_only %>%
# filter(heterosc %in% c("low")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(sex_eff_size = factor(sex_eff_size,
# levels = c("none", "small", "large"))) %>%
ggplot(aes(x = treat_es2, y = ssize, colour = heterosc)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line() +
theme_bw() +
ylab("Sample size") +
xlab("Main treatment effect") +
scale_colour_manual(name = "Effect", values = cbb_palette) + coord_trans(y = "log10") +
theme(text = element_text(size=16))
# plot2
plot2.1 <- all_sex_effs_sex_only %>%
# filter(heterosc %in% c("low")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(sex_eff_size = factor(sex_eff_size,
# levels = c("none", "small", "large"))) %>%
ggplot(aes(x = treat_es2, y = ssize, colour = heterosc)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line() +
theme_bw() +
ylab("Sample size") +
xlab("Main treatment effect") +
scale_colour_manual(name = "Heteroscedasticity", values = cbb_palette) +
xlim(0.3, 0.7) + ylim(5, 25)
# plot2.1Scenario 3: same as no. 1 but with interaction of magnitude present, with the less variable group responding more strongly to the treatment.
alpha <- 0
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -0.5
sex_sd <- 1+alpha
all_effs_sex_large_noh <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
)
)
calc_power_sex_large_noh <- all_effs_sex_large_noh %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(ix_eff_size = rep("none")) %>%
mutate(heterosc = rep("none"))
alpha <- 1
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -0.5
sex_sd <- 1+alpha
all_effs_sex_large_lowh <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
)
)
calc_power_sex_large_lowh <- all_effs_sex_large_lowh %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(ix_eff_size = rep("none")) %>%
mutate(heterosc = rep("low"))
alpha <- 2
no_per_gp <- 5
variable_mean <- 1
variable_sd <- 0.5
treat_es <- 0
sex_es <- 1
female_es <- 0
male_es <- 0
treat_es2 <- treat_es
sex_es2 <- sex_es
ix_es2 <- -0.5
sex_sd <- 1+alpha
all_effs_sex_large_largeh <- rbind(
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.1, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.1, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.2, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.2, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.3, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.3, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.4, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.4, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.5, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.5, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.6, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.6, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.7, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.7, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.8, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.8, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
),
p_value_interaction_crossed(
n_rep = 1000, no_per_gp = no_per_gp, variable_mean = variable_mean,
variable_sd = variable_sd, treat_es = 0.9, sex_es = sex_es,
male_es = male_es, female_es = female_es,
treat_es2 = 0.9, sex_es2 = sex_es2, ix_es2 = ix_es2, sex_sd = sex_sd,
fix_power = T
)
)
calc_power_sex_large_largeh <- all_effs_sex_large_largeh %>%
filter(!Effect == "Residuals ") %>%
group_by(Effect, treat_es2, sex_es2, method, heterosc, ssize, power) %>%
summarise_all(.funs = list(~ sum(. < 0.05, na.rm = TRUE) / n())) %>%
mutate(sex_eff_size = rep("large")) %>%
mutate(ix_eff_size = rep("none")) %>%
mutate(heterosc = rep("large"))
all_sex_effs_sex_only <- rbind(
calc_power_sex_large_noh,
calc_power_sex_large_lowh,
calc_power_sex_large_largeh
)
all_sex_effs_sex_only <- as.data.frame(lapply(all_sex_effs_sex_only, unlist))
# plot1 <- all_sex_effs_sex_only %>%
# filter(Effect %in% c("Treatment")) %>%
# filter(heterosc %in% c("low")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(method = recode(method,
# "ANOVA" = "Factorial")) %>%
# ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
# geom_line() +
# theme_classic() +
# ylab("Statistical power") +
# xlab("Main treatment effect") +
# scale_colour_manual(name = "Sample size", values = cbb_palette) +
# # scale_linetype_discrete(name = "Heteroscedasticity") +
# ylim(0, 1)
# plot2 <- all_sex_effs_sex_only %>%
# filter(Effect %in% c("Treatment")) %>%
# filter(heterosc %in% c("large")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(method = recode(method,
# "ANOVA" = "Factorial")) %>%
# ggplot(aes(x = treat_es2, y = p_value, colour = sample_size)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
# geom_line() +
# theme_classic() +
# ylab("Statistical power") +
# xlab("Main treatment effect") +
# scale_colour_manual(name = "Sample size", values = cbb_palette) +
# # scale_linetype_discrete(name = "Heteroscedasticity") +
# ylim(0, 1)
# ggarrange(plot1, plot2, ncol = 2, nrow = 1, common.legend = TRUE, legend = "right")
plot3 <- all_sex_effs_sex_only %>%
# filter(heterosc %in% c("low")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(sex_eff_size = factor(sex_eff_size,
# levels = c("none", "small", "large"))) %>%
ggplot(aes(x = treat_es2, y = ssize, colour = heterosc)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line() +
theme_bw() +
ylab("Sample size") +
xlab("Main treatment effect") +
scale_colour_manual(name = "Effect", values = cbb_palette) + coord_trans(y = "log10") +
theme(text = element_text(size=16))
# plot3
plot3.1 <- all_sex_effs_sex_only %>%
# filter(heterosc %in% c("low")) %>%
# mutate(sample_size = factor(sample_size,
# levels = c("5M+5F", "10M+10F", "50M+50F"))) %>%
# mutate(sex_eff_size = factor(sex_eff_size,
# levels = c("none", "small", "large"))) %>%
ggplot(aes(x = treat_es2, y = ssize, colour = heterosc)) +
# geom_hline(yintercept = 0.8, linetype = "dashed", colour = "gray") +
geom_line() +
theme_bw() +
ylab("Sample size") +
xlab("Main treatment effect") +
scale_colour_manual(name = "Heteroscedasticity", values = cbb_palette) +
xlim(0.3, 0.7) + ylim(5, 2200)
# plot3.1(In all below examples the y axis (sample size) indicates the number of individuals per group and sex - hence the total sample size is 4 times the number indicated on the y axis.)
Combined plot for scenarios of no interaction (A), interaction with stronger effect in more variable sex (B) and stronger effect in less variable sex (C) for full range of effect sizes and on log-scale for sample size:
comb1 <- ggarrange(plot1, plot2, plot3, ncol = 1, nrow = 3, common.legend = TRUE, legend = "bottom", labels = c("A", "B", "C"))
comb1Combined plot for scenarios of no interaction (A), interaction with stronger effect in more variable sex (B) and stronger effect in less variable sex (C) for moderate-sized effect sizes using identity scale for sample size:
comb2 <- ggarrange(plot1.1, plot2.1, plot3.1, ncol = 1, nrow = 3, common.legend = TRUE, legend = "bottom", labels = c("A", "B", "C"))
comb2Combined plots demonstrating the impact of heteroscedasticity on power (repeated from section Section 1.2) for a no-interaction case (A) and interaction-of-magnitude case (B; stronger effect in the more variable sex):
comb3 <- ggarrange(plotA, plotB, ncol = 2, nrow = 1, common.legend = TRUE, legend = "bottom", labels = c("A", "B"))
comb3Below plots show hypothetical scenarios when heterogeneity is a natural outcome of the underlying statistical process. Diagonal lines connect averages belonging to one treatment group (control vs. experimental). Thick vertical lines show SD of each subgroup.
Simulation draws from \(Pois(\lambda_i)\) assumes sex/treatment specific averages (\(\lambda_i\)) - see lambda1...lambda4 below.
N <- 500
lambda1 <- 3
lambda2 <- 9
lambda3 <- 8
lambda4 <- 15
toy_poisson <- data.frame(
Treatment = c(rep("Ctrl", 2*N), rep("Exp", 2*N)),
Sex = c(rep("F", N), rep("M", N), rep("F", N), rep("M", N)),
Response = c(rpois(N, lambda1), rpois(N, lambda2), rpois(N, lambda3), rpois(N, lambda4))
)
toy_poisson$Group <- interaction(toy_poisson$Sex, toy_poisson$Treatment)
plot1 <- ggplot(data = toy_poisson, aes(x = Response, y = Group)) +
geom_density_ridges(aes(fill = Sex, scale = 0.8),
stat = "binline", bins = 15, alpha = 0.6, draw_baseline = FALSE
) +
theme_ridges() +
geom_point(size = 4, stat = "summary", fun = "mean") +
stat_summary(
fun = "mean", geom = "line", linewidth = 1,
aes(group = Treatment, linetype = Treatment)
) +
stat_summary(
fun.max = function(x) mean(x) + sd(x),
fun.min = function(x) mean(x) - sd(x), geom = "errorbar",
linewidth = 1.5, width = 0
) +
theme_bw() +
coord_flip() +
theme(text = element_text(size = 16))
# plot1Simulation samples from \(Lognormal(\mu_i, \sigma^2=0.1)\) with sex/treatment specific means (mu1...mu4) below.
N <- 500
mu1 <- 0.01
mu2 <- 0.5
mu3 <- 0.4
mu4 <- 0.8
sd <- 0.31
toy_lognorm <- data.frame(
Treatment = c(rep("Ctrl", 2*N), rep("Exp", 2*N)),
Sex = c(rep("F", N), rep("M", N), rep("F", N), rep("M", N)),
Response = c(rlnorm(N, mu1, sd), rlnorm(N, mu2, sd), rlnorm(N, mu3, sd), rlnorm(N, mu4, sd))
)
toy_lognorm$Group <- interaction(toy_poisson$Sex, toy_poisson$Treatment)
plot2 <- ggplot(data = toy_lognorm, aes(x = Response, y = Group)) +
geom_density_ridges(aes(fill = Sex, scale = 0.8),
stat = "binline", bins = 15, alpha = 0.6, draw_baseline = FALSE
) +
theme_ridges() +
geom_point(size = 4, stat = "summary", fun = "mean") +
stat_summary(
fun = "mean", geom = "line", linewidth = 1,
aes(group = Treatment, linetype = Treatment)
) +
stat_summary(
fun.max = function(x) mean(x) + sd(x),
fun.min = function(x) mean(x) - sd(x), geom = "errorbar",
linewidth = 1.5, width = 0
) +
theme_bw() +
coord_flip() +
theme(text = element_text(size = 16))
# plot2ggarrange(plot1, plot2, ncol = 2, nrow = 1,
common.legend = TRUE, legend = "bottom", labels = c("A", "B"))Heterogeneity in variances may also arise if the underlying phenotypic mapping function is non linear and has saturating/desaturating character in a range modified by sex dimorphism and/or experimental treatment.
Here we assume a logistic function:
\[f(x)=\frac{L}{e^{-k(x-x_0)}}\]
with parameter values: L = 3.5, k = 1, x0 = 1.
mapping_logistic <- function(x, L = 3.5, x0 = 1, k = 1) L/(1+exp(-1*k*(x - x0)))Using this function, we can generate data assuming some treatment and sex differences:
SD_N <- 0.5
mu_env1 <- 1
mu_env2 <- 2
mu_env3 <- 3
mu_env4 <- 4.5
toy_nonlin <- data.frame(
Treatment = c(rep("Ctrl", 2 * N), rep("Exp", 2 * N)),
Sex = c(rep("F", N), rep("M", N), rep("F", N), rep("M", N)),
Env = c(
rnorm(N, mu_env1, SD_N), rnorm(N, mu_env2, SD_N),
rnorm(N, mu_env3, SD_N), rnorm(N, mu_env4, SD_N)
)
)
toy_nonlin$Response <- mapping_logistic(toy_nonlin$Env)ggplot() +
geom_point(
data = toy_nonlin,
aes(y = Env, x = Response, colour = Sex),
pch = 1, size = 3, alpha = 0.8
) +
geom_line(
colour = "black", linetype = 2,
aes(y = seq(-1, 8, 0.1), x = mapping_logistic(seq(-1, 8, 0.1)))
) +
geom_density_ridges(
data = toy_nonlin %>%
group_by(Sex, Treatment) %>%
mutate(meanEnv = mean(Env)),
aes(x = Response, group = meanEnv, y = meanEnv, fill = Sex),
stat = "binline", alpha = 0.8, draw_baseline = F
) +
geom_point(
data = toy_nonlin %>%
group_by(Sex, Treatment) %>%
summarise(meanEnv = mean(Env), meanResp = mean(Response)),
aes(x = meanResp, y = meanEnv),
size = 4
) +
geom_line(
data = toy_nonlin %>%
group_by(Sex, Treatment) %>%
summarise(meanEnv = mean(Env), meanResp = mean(Response)),
aes(x = meanResp, y = meanEnv, linetype = Treatment),
linewidth = 1
) +
geom_errorbarh(
data = toy_nonlin %>%
group_by(Sex, Treatment) %>%
summarise(meanEnv = mean(Env), meanResp = mean(Response),
minResp = mean(Response) - sd(Response),
maxResp = mean(Response) + sd(Response)),
aes(y = meanEnv, xmin = minResp, xmax = maxResp),
linewidth = 1.5, height = 0
) +
# stat_summary(data = toy_nonlin %>%
# group_by(Sex, Treatment) %>%
# mutate(meanEnv = mean(Env)),
# fun = "mean", geom = "line", linewidth = 1,
# aes(x = Response, group = meanEnv, y = meanEnv, linetype = Treatment)
# ) +
coord_flip() + ylab("Manipulated environmental variable") +
theme_ridges() +
theme_bw() + theme(text = element_text(size = 16))Here we use previously published synthetic data to demonstrate the degree and ubiquity of mean-variance relationships typical to the Taylor’s law.
First let’s load the relevant data.
data_all <- read_excel(here("R", "ALL_TRAITS_20200514.xlsx"), sheet = "main")
# names(data_all)
## analyses on adiposity [g] data only:
data_all %>%
filter(Trait == "Adiposity") %>%
count(Paper_ID)# A tibble: 12 × 2
Paper_ID n
<chr> <int>
1 Adedeji2019 4
2 Armitage2007 2
3 Barbosa2020 2
4 CastroBarbosa2016 16
5 Ding2014 1
6 King2013 21
7 Li2012 1
8 MartinsTerra2019 4
9 Masuyama2015 2
10 Sarker2018 12
11 Tait2015 4
12 Winther2019 2
# offspring ages
summary(data_all$Age_at_Measurement_Days[data_all$Trait == "Adiposity"]) Min. 1st Qu. Median Mean 3rd Qu. Max.
84.0 84.0 91.0 126.8 180.0 270.0
# species
summary(factor(data_all$Rodent_Type[data_all$Trait == "Adiposity"]))Mouse Rat
37 34
# merge control and treatment data (all traits)
dataA <- tibble(
mean = c(data_all$Mean_Treatment, data_all$Mean_Control),
sd = c(data_all$SD_Treatment, data_all$SD_Control),
size = c(data_all$Sample_Size_n_Treatment, data_all$Sample_Size_n_Control),
sex = c(data_all$Offspring_Sex, data_all$Offspring_Sex),
trait = c(data_all$Trait, data_all$Trait),
species = c(data_all$Rodent_Type, data_all$Rodent_Type)
)Plotting individual subplots for mice, rats, and both species jointly.
dataA %>%
filter(trait == "Adiposity" & species == "Mouse") %$%
cor.test(log(mean), log(sd)) -> cor_mice # calculate correlation
figA <- dataA %>%
filter(trait == "Adiposity") %>%
filter(species == "Mouse") %>%
ggplot(aes(x = log(mean), y = log(sd), col = sex)) +
geom_point(aes(size = size), alpha = 0.5) +
geom_smooth(method = "lm", se = FALSE) +
# coord_trans(x = "log10", y = "log10") +
theme_bw() +
theme(legend.position = "none") +
labs(title = "A. Mouse adiposity", x = "ln (mean)", y = "ln (SD)") +
geom_text(x = -5.2, y = 1.5, aes(label = sprintf("italic(r) == %.2f", cor_mice$estimate)), parse = TRUE, hjust = 0, col = "black", size = 5)
dataA %>%
filter(trait == "Adiposity" & species == "Rat") %$%
cor.test(log(mean), log(sd)) -> cor_rats # calculate correlation
figB <- dataA %>%
filter(trait == "Adiposity") %>%
filter(species == "Rat") %>%
ggplot(aes(x = log(mean), y = log(sd), col = sex)) +
geom_point(aes(size = size), alpha = 0.5) +
geom_smooth(method = "lm", se = FALSE) +
# coord_trans(x = "log10", y = "log10") +
theme_bw() +
theme(legend.position = "none") +
labs(title = "B. Rat adiposity", x = "ln (mean)", y = "ln (SD)") +
geom_text(x = -1, y = 4, aes(label = sprintf("italic(r) == %.2f", cor_rats$estimate)), parse = TRUE, hjust = 0, col = "black", size = 5)
# correlation - all
dataA %>%
filter(trait == "Adiposity") %$%
cor.test(log(mean), log(sd)) -> cor_all
figC <- dataA %>%
filter(trait == "Adiposity") %>%
ggplot(aes(x = log(mean), y = log(sd), col = sex)) +
geom_point(aes(size = size), alpha = 0.5) +
geom_smooth(method = "lm", se = FALSE) +
# coord_trans(x = "log10", y = "log10") +
theme_bw() +
theme(legend.position = "bottom") +
labs(title = "C. Mouse and rat adiposity", x = "ln (mean)", y = "ln (SD)") +
geom_text(x = -5, y = 4, aes(label = sprintf("italic(r) == %.2f", cor_all$estimate)), parse = TRUE, hjust = 0, col = "black", size = 5)
# assemble the panels using patchwork package
figure1 <- figA / figB / figC +
plot_layout(ncol = 1, nrow = 3, heights = c(1, 1, 1.1)) #+
# plot_annotation(tag_levels = "A")
# ggsave(plot = figure1, "Figure1_v0.png", width = 8, height = 14, units = "cm", dpi = "retina", scale = 2)
# ggsave(plot = figure1, "Figure1_v0.pdf", width = 8, height = 14, units = "cm", scale = 2)
figure1